Two approaches

1. The latent variable interpretation. Somewhat abstract, but simple enough.

2. The odds-ratio and approximate relative risk approach. Only for logistic regression though!

The odds

• Let \(p = P(Y=1)\). Recall that \(P(Y=1) = F(\beta^T X)\) in a binary regression model.
• The odds for an event with probability \(p\) is defined as \(\text{odds} = p/(1-p)\).
• In logistic regression, \(p = \frac{\exp[\beta^{T}x]}{1+\exp[\beta^{T}x]}\).
• Hence the odds are \[\frac{\frac{\exp[\beta^{T}x]}{1+\exp[\beta^{T}x]}}{1-\frac{\exp[\beta^{T}x]}{1+\exp[\beta^{T}x]}} = \frac{\frac{\exp[\beta^{T}x]}{1+\exp[\beta^{T}x]}}{\frac{1}{1+\exp[\beta^{T}x]}} = \exp[\beta^Tx].\]

The odds-ratio

• Let \(p\) and \(q\) be two probabilities. The odds-ratio is the ratio of odds with probabilities \(p\) and \(q\). \[\text{OR} = \frac{p/(1-p)}{q/(1-q)}.\]
• In the logistic regression model, the odds-ratio between two covariates \(x'\) (for \(p\)) and \(x\) (for \(q\)) is \[\frac{p/(1-p)}{q/(1-q)} = \frac{\exp[\beta^Tx']}{\exp[\beta^Tx]} = \exp{[\beta^T(x'-x)]}.\]
• Now suppose that \(x' = x\) everywhere except on the \(j\)th coordinate, where it equals \(x+1\). Then \(\beta^T(x'-x) = \beta_j\).

The log-odds and relative risk

• From the equation for the odds-ratio we find that \[\frac{p/(1-p)}{q/(1-q)} = \exp{[\beta^T(x'-x)] = \exp{[-\beta_j]}}.\]
• Hence the logarithm of the odds-ratio equals \(\beta_j\)!
• The relative risk of two events is \(\text{RR} = p/q\). How much more probable is \(Y\) under \(x'\) than \(Y\) under \(x\)? Much easier to interpret than the odds ratio.
• When \(p,q\) are small, the odds ratio approximates the relative risk.
• So when the probabilities are small, \(\exp{[\beta_j]\approx \text{RR}}\) of increasing the value of \(x_j\) by \(1\).

Example

mod = smf.probit("admit ~ gre + gpa + rank", data=admission).fit()

Optimization terminated successfully.
         Current function value: 0.574351
         Iterations 5

To intepret the coefficients we must exponentiate them.

np.exp(mod.params)

Intercept    0.123501
gre          1.001399
gpa          1.590995
rank         0.717694
dtype: float64

Hence increasing gpa with 1 increases the probability of being admitted by \(60\%\), everything else held equal. Increasing the rank of the university by one decreases the probability by \(30\%\) though.

The multivariate normal

\[f(\boldsymbol{x};\boldsymbol{\mu},\Sigma)=\frac{\exp\left(-\frac{1}{2}\left(\boldsymbol{x}-\boldsymbol{\mu}\right)^{T}\Sigma^{-1}\left(\boldsymbol{x}-\boldsymbol{\mu}\right)\right)}{\sqrt{(2\pi)^{k}|\Sigma|}}\]

• \(\boldsymbol{x}\in\mathbb{R}^k\) is the vector of observations

• \(\Sigma\) is the \(k\times k\) covariance matrix.

• \(\boldsymbol{\mu}\in\mathbb{R}^k\) is the vector of means. * \(|A|\) denotes the determinant of \(A\).

The multivariate central limit theorem

• Recall that \(\sqrt{n}(\overline{x}-\mu) \to N(0, \sigma^2)\), where \(\sigma^2 = \operatorname{Var}(X)\) when \(x\) is univariate and the mean is obtained from iid observations.
• When \(\boldsymbol{X_1}, \boldsymbol{X_2},\ldots, \boldsymbol{X_k}\) is a sequence of multivariate observations, define its multivariate mean \(\boldsymbol{X}\) component-wise, i.e., \(\overline{X_{\cdot 1}}, \overline{X_{\cdot 2}}, \ldots, \overline{X_{\cdot k}}\).
• Remember that the covariance matrix of a random variable is \(E[(X-\mu)(X-\mu)^T]\). It is a positive definite matrix that summarizes the dependence between the variables.

An example

• Let \(U\), \(V\), \(W\) be independent uniforms. These have means \(1/2\) and variances \(1/12\).
• Define \(X_1 = U+V\) and \(X_2 = U+V+W\).
• Then the (co)variances are:
  • \(\operatorname{Var}(X_1) = 1/6\)
  • \(\operatorname{Var}(X_2) = 1/4\)
  • \(\operatorname{Cov}(X_1, X_2) = 1/6\)
• Question: What are the means of \(X_1\) and \(X_2\)?
  • Answer: \(1\) and \(3/2\).
• Question: What is the covariance matrix of \(\sqrt{n}(\overline{X_1}, \overline{X_2})\) when \(n\to\infty\)?
  • Answer: The covariance matrix of \((X_1, X_2)\), by the central limit theorem.

The Fisher information matrix

• The Fisher information is \[I(\theta)=E\left[\left(\frac{\partial}{\partial\theta}\log f(X;\theta\right)^{2}\mid\theta\right]\]
• Almost always equal to \[I(\theta)=-E\left[\frac{\partial^{2}}{\partial\theta^{2}}\log f(X;\theta)\mid\theta\right]\]
• Summarizes the information about the parameter contaied in one observation.

Example Fisher information

• Exponential distribution: \(f(x;\theta)=\theta e^{-\theta x}\).
• \(\log f(x;\theta)=\log\theta-\theta x\)
• \(\frac{\partial}{\partial\theta}\log f(x;\theta)=\frac{1}{\theta}-x\)
• \(\frac{\partial^{2}}{\partial\theta^{2}}\log f(x;\theta)=-\frac{1}{\theta^{2}}\)
• Hence \(I(\theta)=-E\left[\frac{\partial^{2}}{\partial\theta^{2}}\log f(X;\theta)\right]=\frac{1}{\theta^{2}}\).

The asymptotic distribution of the maximum likelihood estimator

• If \(\theta\) is a vector, the multivariate Fisher information is the matrix with elements \[I(\theta)_{ij} = -E\left[\frac{\partial^{2}}{\partial\theta_i\theta_j}\log f(X;\theta)\mid\theta\right]\]

• Question: Why do we care about this?
  • Answer: You can construct confidence intervals, hypothesis tests, and you learn how quickly the parameter estimates approach their true values.

Finding the asymptotic covariance matrix in Python

The estimated covariance matrix of the parameter estimates an be found by calling

mod.cov_params()

• This is approximately equal to \(\frac{1}{n}{I^{-1}(\theta)}\), where \(\theta\) is the true parameter value.
• Don’t misinterpret this! It’s the covariance of the parameter estimates, not the data itself!

An example contingency table

Model: \(P(\text{Urinary tract infection}) = F(\beta_0+\beta_11[\text{Diaphragm use}]),\) where \(F\) is equal to, e.g., the logistic CDF.

• What is \(\hat{\beta}_0\)?
  • Since the maximum likelihood estimator of \(\hat{p}\) in a Bernoulli model is the mean, it is \(140/(140 + 290) \approx 0.33\).
• What is \(\hat{\beta}_1\)?
  • This part only invoves “Yes” events! We want to maximize the likelihood, so it’s best to make the probability equal to \(1\), which implies \(\beta_1 = \infty\)!

Perfect separation

Example of perfect separation

Definition

• Let \(X\) be the model matrix, i.e., the matrix of covariates and the intercept.
• Define \(X_0\) as the submatrix where the outcome is \(Y=0\) and \(X_1\) the submatrix where \(Y=1\).

More detailed explanations can be found here.

Example data (i)

A researcher wishes to know if caffeine improves performance on a memory test. Volunteers consume different amounts of caffeine from 0 to 500 mg, and their score on the memory test is recorded. Source.

Example data (ii)

• This data is on a different form than usual! We have n giving us the number of observations with the given covariate levels and a the numbe who got an A!
• This can be called a binomial regression model. It’s equivalent to binary regression models, but with a more compact representation.

Fitting the model

import statsmodels.api as sm
model = smf.glm("a + I(n - a) ~ caffeine", family=sm.families.Binomial(), data=data).fit()
model.params

Intercept    0.238469
caffeine    -0.006442
dtype: float64

Is the model good?

import matplotlib.pyplot as plt
logistic = lambda x: np.log(x) - np.log(1-x)
plt.plot(data.caffeine, logistic(model.predict()))
plt.scatter(data.caffeine, logistic(data.prob))

<matplotlib.collections.PathCollection at 0x1f8a99ac490>

• Question: What’s going on here?
  • There seems to be a \(U\)-shape! Some caffeine is good, but more is bad.

How well does it fit?

plt.scatter(model.predict(), data.prob)
plt.plot([0,0.6], [0, 0.6])

• These points should, ideally, match the straight line!
• But is the discrepancy just due to uncertainty?

Trying another model: Adding a log term.

model_log = smf.glm(
  "a + I(n - a) ~ caffeine + np.log(caffeine + 1)", 
  family=sm.families.Binomial(), 
  data=data).fit()
plt.scatter(model_log.predict(), data.prob, color = "blue")
plt.scatter(model.predict(), data.prob, color = "red")
plt.plot([0,0.6], [0, 0.6])

What does it predict?

<matplotlib.collections.PathCollection at 0x1f8aa144a30>

• Question: It looks like the model with the log term fits the data better! But what might go wrong?
  • Answer: Maybe it’s just noise! Maybe the low-caffeine observations are too low, or the the medium-caffeine observations too high.

The AIC: Definition

• Akaike’s information criterion (AIC), defined as \[\text{AIC} = 2p - 2l,\] where \(p\) is the number of estimated parameters and \(l\) is the log-likelihood at the estimated parameters.
• The smaller the AIC the better!
• We add \(2p\) since a larger number of parameters allows for more noise. It penalizes model complexity.
• Penalization of model complexity is always needed when comparing model with different number of parameters!
• There is a theoretical argument for the choice \(p\)! but that is beyond the scope of this course.

The AIC: Application

model.aic

55.87004439904055

model_log.aic

45.844885174715465

• Hence the logarithm term is worth it. This is a large reduction in AIC; you’ll get an intuition about this after a while.
• Key take-away: Drinking some coffee will help your grades!

The AIC: Why and when?

• Question: Why use the AIC?
  • It is simple, it is well-known, and it performs well.
• You can compare AIC across different models!

model_log_probit = smf.glm(
  "a + I(n - a) ~ caffeine + np.log(caffeine + 1)", 
  family=sm.families.Binomial(sm.genmod.families.links.probit()), 
  data=data).fit()
model_log_cauchit = smf.glm(
  "a + I(n - a) ~ caffeine + np.log(caffeine + 1)", 
  family=sm.families.Binomial(sm.genmod.families.links.cauchy()), 
  data=data).fit()
model_log_cloglog = smf.glm(
  "a + I(n - a) ~ caffeine + np.log(caffeine + 1)", 
  family=sm.families.Binomial(sm.genmod.families.links.cloglog()), 
  data=data).fit()  
{"logit" : model_log.aic, "probit" : model_log_probit.aic, "cauchit" : model_log_cauchit.aic, "cloglog":  model_log_cloglog.aic}

{'logit': 45.844885174715465,
 'probit': 45.37382082095568,
 'cauchit': 52.37445834568989,
 'cloglog': 46.901711989146094}

\(0-1\) predictions

• The predict method gives you a probability.
• But often you have to make a make a prediction of \(\hat{Y}\) itself, e.g., will you follow a marketing lead? Do you believe candidate will get an \(A\) when he drinks 5 cups of coffee?
• You usually construct predictions of \(Y\) using a threshold: \[\hat{Y}=1 \iff \hat{p} > c,\] where \(\hat{p}\) is a predicted probability obtained using predict and \(\hat{Y}\) is the predicted value of \(Y\).
• For instance, you could decide to do the action if if \(\hat{p} > 0.5\).

Admission data (again)

admission = pd.read_csv("https://stats.idre.ucla.edu/stat/data/binary.csv")
admission.head()
model = smf.logit("admit ~ gre + gpa + C(rank)", data = admission).fit()
plt.hist(model.predict())

Optimization terminated successfully.
         Current function value: 0.573147
         Iterations 6

(array([31., 69., 62., 67., 66., 40., 26., 20.,  8., 11.]),
 array([0.05878643, 0.12674861, 0.19471079, 0.26267297, 0.33063516,
        0.39859734, 0.46655952, 0.5345217 , 0.60248388, 0.67044606,
        0.73840825]),
 <BarContainer object of 10 artists>)

Confusion matrix

The confusion matrix shows shows true positives, false negatives, false positives, and true negatives.

Confusion matrix for admission with c = 0.5

model.pred_table(threshold = 0.5)

array([[254.,  19.],
       [ 97.,  30.]])

• Question: How many true positives are there?
  • Answer: True positives happen when both \(\hat{Y}=1\) and \(Y=1\). Hence the correct answer is \(254\).
• Question: How many false negatives are there?
  • Answer: False negatives correspond to \(\hat{Y}=0\) and \(Y=1\). There are 19 of these.

Confusion matrix for admission with c = 0.6

model.pred_table(threshold = 0.6)

array([[266.,   7.],
       [114.,  13.]])

model.pred_table(threshold = 0.5)

array([[254.,  19.],
       [ 97.,  30.]])

• Question: How can you make sure the number of false positives is \(0\)?
  • Answer: Make the threshold \(c=1\)!

True positive rate and false positive rate

\[ \text{true positive rate}=\frac{\text{true positives}}{\text{true positives} + \text{false negatives}} \] Also known as sensitivity, recall, and hit rate. Abbreviated as tpr.

\[ \text{false positive rate}=\frac{\text{false positives}}{\text{false positives} + \text{true negatives}} \]

Also known as fall-out. Abbreviated as fpr.

• Here \(TPR = \frac{30}{30 + 97} = 0.24\) and \(FPR = \frac{19}{19 + 254} = 0.07\).

Receiver operating characteristic curve

The receiver operating characteristic curve. Source: Wikipedia.

How to make these plots

from sklearn import metrics
fpr, tpr, _ = metrics.roc_curve(admission.admit, model.predict())
plt.plot(fpr, tpr)
plt.show()

The area under the curve (i)

The area under the curve (ii)

• The area under the curve (AUC) is the integral of the receiver operating characteristic curve.
  • \(0.5 \leq AUC \leq 1\).
  • \(AUC = 1\) if we can predict perfectly.
  • \(AUC = 0.5\) if we can’t predict at all.
• Extremely widespread!
• Since the AUC doesn’t penalize model complexity it shouldn’t be used to compare models!

The \(R^2\) for binary regression models

Remember the setup for the \(R^{2}\) for non-linear regression: \[ 1-\frac{\sum_{i=1}^{n}(y_{i}-f(x;a_{1},\ldots a_{p}))^{2}}{\sum_{i=1}^{n}(y_{i}-m)^{2}}\] In general we could write \[R^{2}=1-\frac{R(x)}{R(m)}\] Where

• \(R(m)\) measures the distance between our best predictions and the data \(y\), and
• \(R(x)\) measures the distance between our best predictions and the data \(y\) when we take \(x\) into account.

The log-likelihood and the \(R^2\)

We use the log-likelihood instead of the quadratic loss! So, letting \(p_i = f(x;a_{1},\ldots a_{p})\), define the fitted likelihood \[l(x) = R(x) = \sum_{i=1}^{n}-y_{i}\log \hat{p}_{i}-(1-y_{i})\log(1-\hat{p}_{i}).\]

\[l(m) = R(m) = \sum_{i=1}^{n}-y_{i}\log m-(1-y_{i})\log(1-m),\] is the log likelihood of the null model. It measures how well we can predict \(y\) when we don’t know any \(x\) at all.

McFadden’s \(R^2\)

Now define the pseudo-\(R^2\) or McFadden \(R^2\) as \[R^2_{\textrm{McFadden}} = 1 - \frac{l(x)}{l(m)}.\]

• Since McFadden’s \(R^2\) doesn’t penalize model complexity it shouldn’t be used to compare models!

An example

probit_fit = smf.probit("admit ~ gre + C(rank) + gpa", data = admission).fit()
logit_fit = smf.logit("admit ~ gre + C(rank) + gpa", data = admission).fit()

def rsq_mcfadden(fit):
    lower = fit.llnull
    upper = fit.llf
    return 1 - upper / lower

print(rsq_mcfadden(probit_fit), rsq_mcfadden(logit_fit))

Optimization terminated successfully.
         Current function value: 0.573016
         Iterations 5
Optimization terminated successfully.
         Current function value: 0.573147
         Iterations 6
0.08313059666890721 0.08292194470084713

• Thus the \(R^2\)s are \(0.0831\) and \(0.0829\). But we may also access McFadden’s \(R^2\) directly: print(probit_fit.prsquared, logit_fit.prsquared).

• Question: Can you compare these \(R^2\)s to the \(R^2\)s from linear regression?

  • Answer: No, as we are looking the likelihoods now, not the quadratic distance. Also, the “ordinary” \(R^2\) tends to be small for binary data.

McFadden’s \(R^2\) vs the ordinary \(R^2\)

The following picture, Figure 5 from Chapter 5 of McFadden’s Urban Travel Demand: A Behavioral Analysis (1996) illustrates the typical relationsship between the least squares \(R^2\) (on the \(y\)-axis) and the pseudo-\(R^2\) on the \(x\)-axis.

McFadden’s comparison if the ordinary R^2 and McFadden’s R^2. Source: McFadden, D. Urban Travel Demand: A Behavioral Analysis (1996)

Summary part 1

1. The parameters of a logistic regression model can be interpreted as the log-odds-ratio. The log-odds-ratio approximates the relative risk.
2. Under general conditions, the asymptotic distribution of a maximum likelihood estimator \(\hat{\theta}\) is \(\sqrt{n}\left(\hat{\theta}-\theta\right)\to N(0,I(\theta)^{-1})\), where \(I(\theta)\) is the Fisher information matrix.
3. We can use this to calculate confidence intervals and construct hypothesis tests for the parameters of a binary regression model.
4. The maximum likelihood estimates \(\hat{\theta}\) of a binary regression model are finite if and only if the data isn’t quasi-perfectly separated.

Summary part 2

5. Akaike’s information criterion (AIC) can be used to compare models with different covariates or link functions.
6. To predict \(Y\) from fitted probabilities you need to choose a threshold. Each threshold is associated with a confusion matrix, displaying the relationsship between \(\hat{Y}\) and \(Y\).
7. The ROC displays the relationsship between the true positive rate and the false negative rate.
8. The AUC measures how good your model is at predicting, with \(0.5\) being no better than chance and \(1\) being perfect. The AUC should not be used to compare models. Use the AIC for that.
9. McFadden’s \(R^2\) is a generalization of the \(R^2\) to binary models.